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ABSTRACT 


The  validity  of  linear  smoothing  of  topography  for  numerical  weather 
prediction  and  the  variation  of  mountain  drag  with  mountain  height  and 
static  stability  are  examined  in  this  study.   In  the  model  a  constant 
geostrophic  current  is  perpendicular  to  the  mountain  range  and  the 
height  of  the  mountain  is  independent  of  y.   The  hydrostatic  Boussinesq 
equations  are  used  with  motion  bounded  at  the  top  by  a  rigid  plane  at 
z  =  D.   A  modified  coordinate  system  similar  to  Phillips'  sigma  system 
was  used.   Solutions  were  obtained  using  a  smoothed  mountain  profile. 
These  solutions  were  compared  with  smoothed  solutions  obtained  from  the 
unsmoothed  mountain.   The  comparison  of  these  solutions  shows  that  an 
error  is  introduced  when  non- linear  terms  become  sufficiently  large. 
Values  of  the  mountain  drag  for  differing  values  of  mountain  height  at 
a  given  static  stability  and  for  differing  values  of  static  stability 
at  a  given  mountain  height  were  computed.   Mountain  drag  was  found  to 
vary  quadratically  with  mountain  height  and  linearly  with  static 
stability. 
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I.   INTRODUCTION 

The  use  of  large  scale  numerical  forecast  models  makes  it  necessary 
to  introduce  a  topographical  model  of  the  surface  of  the  earth  for  use 
in  the  lower  boundary  condition.   Cressman  [1960]  notes  that  since  the 
topography  of  the  surface  of  the  earth  has  many  small  scale  features, 
it  is  necessary  to  smooth  the  topography  when  it  is  used  as  a  boundary 
condition  in  a  numerical  forecast  model.   Mintz  [1956]  and  Cressman 
[1960]  indicate  that  the  influence  of  small  scale  features  of  the  topog- 
raphy can  be  included  in  a  surface  drag  coefficient  which  varies  from 
point  to  point. 

One  purpose  of  this  study  is  to  examine  the  validity  of  linear 
smoothing  of  the  topography  in  numerical  weather  prediction.   To  accom- 
plish this  the  flow  over  a  small  scale  mountain  was  computed  using  a 
small  grid  mesh.   The  results  of  these  computations  were  then  smoothed. 
The  mountain  was  then  smoothed  and  the  solution  was  recomputed.   Since 
the  parameters  predicted  in  large-scale  forecasts  are  interpreted  as 
area  averages,  a  comparison  of  the  two  results  should  show  little,  if 
any,  difference  if  the  smoothing  technique  is  valid. 

Mountain  drag,  or  mountain  torque,  which  arises  from  a  difference 

in  pressure  on  the  two  sides  of  a  mountain,  is  important  in  the  general 

circulation  of  the  atmosphere.   Studies  of  mountain  drag  associated  with 

smaller  scale  mountains  could  be  used  to  improve  the  calculation  of  the 

surface  drag  coefficient  which  was  mentioned  above.   Miles  [1968]  showed 

that  wave  drag  on  an  obstacle  is  a  function  of  the  height  of  the  obstacle 

and  of  the  mean  flow  perpendicular  to  the  obstacle.   The  present  paper 

investigates  the  effect  of  mountain  height  and  initial  static  stability 

on  mountain  drag. 
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This  study  employes  a  two-dimensional,  primitive  equation  model. 
The  geostrophic  current  perpendicular  to  the  mountain  is  constant. 
The  mountain  range  is  infinite  in  lateral  extent  which  makes  the  model 
independent  of  y.   A  special  coordinate  system  similar  to  that  developed 
by  Phillips  [1957]  was  used.   The  equations  were  solved  numerically 
with  the  use  of  a  special  finite  difference  scheme  which  was  developed 
by  A.  Arakawa  (see  Langlois  and  Kwok,  [1968]). 
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II.   THE  BASIC  EQUATIONS 

The  basic  equations  used  in  this  study  are  the  Boussinesq  equations, 
These  equations  will  be  applied  to  a  model  in  which  the  flow  is  bounded 
above  by  a  rigid  horizontal  plane.   These  equations  neglect  the  compres- 
sibility of  the  atmosphere  which,  for  this  study,  should  not  be  of 
qualitative  importance.   It  will  be  seen  that  the  internal  scale  of  the 
flow  is  small  compared  with  the  depth  of  the  fluid.   The  hydrostatic 
approximation  is  made  and  heating  and  friction  terms  are  neglected. 
The  Boussinesq  equations  (Ogura  and  Phillips  [1962])  may  be  expressed 
in  the  z-coordinate  system  as 

—  +  V-W  +  w  g  =   -  VI    -    fk  x  V,  (1) 

ft  +  v-™  +  -  i  =  °>  <2> 

o 

7i   .    dw 


V-V  +  ^  =  0  (4) 


where 


e  =  T(po/P)H  -  eo 

In  this  model  the  only  quantity  which  varies  with  latitude  is  §  and  it 
varies  linearly.   The  Coriolis  parameter  is  constant. 

To  simplify  the  finite  differencing  of  the  basic  equations  a 
modified  coordinate  system  is  introduced  in  which  the  lower  boundary 
surface  becomes  a  coordinate.   This  system  makes  it  possible  to  intro- 
duce the  effects  of  topography  without  employing  uncentered  horizontal 
space  differences  in  the  vicinity  of  the  mountains.   It  also  allows 
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application  of  the  lower  boundary  condition  at  the  actual  surface  of 
the  mountain.   In  the  proposed  coordinate  system  the  vertical  advection 
terms  are  identically  zero  at  both  the  upper  and  lower  boundaries.   This 
system  of  coordinates  is  very  similar  to  the  sigma  system  introduced  by 
Phillips  [1957]. 

The  vertical  coordinate,  zeta,  is  defined  as: 

z  -  z 

c  -  -r^  (5) 

where  : 

Z  =  D  -  z 

s 

According  to  this  definition  Q    =  0  at  the  lower  boundary  and  Q   =  1  at 
the  upper  boundary.   Fig.  1  illustrates  this  coordinate  system. 

Transformation  of  equations  (1)  through  (4)  from  the  z-coordinate 
system  into  the  ^-coordinate  system  yields: 


g  +  v.w  +  cf=-™-!(4->h*s-f£xv 


(6) 


|f  +  V-V9  +  C  |l  -  0  (7) 

|  =  f-ze  (8) 


V-ZV  +  Zg=0  (9) 

In  equations  (6)  through  (9)  V  is  the  horizontal  gradient  operator  on 

zeta  surfaces. 

Since  some  of  the  following  computations  involve  vertical  averaging 

it  is  convenient  to  express  the  equations  in  flux  form.   Equations  (6) 

through  (9)  in  component  form  become : 

d(M  )  -w  n 

-gS- +  £<»)=-.  §  +|_H(C-l)f  +fMy         (10) 
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-*-  +  £(v)  =  -  Z  §f  -  fM  (11) 


ht     t  ±y,   -  -  *  ay 


X 


If  +  f '(e)  -  o  (12) 


ax    +  z  ac     ° 


(13) 


(14) 


where  : 


£'(s)  -  |j  (sHx)  +  |j  (C.X) 

S  =  U,  V,  0 

and  M  =  uZ,  M  =  vZ ,  H  =  9Z . 

x      y 

The  boundary  conditions  are: 

C  =  0  at  C  =  0,  and  (15) 

C  =  0  at  C  =  1.  (16) 

For  this  study  the  flow  will  be  independent  of  y  and  periodic  in  x. 

Since  $  is  the  only  variable  allowed  to  change  in  the  y  direction 

■57—  =  0  for  all  other  variables, 
dy 

To  obtain  an  expression  for  $  equation  (13)  may  be  integrated  with 
respect  to  Q  , 

c 

$  =  [I-  H  dC  +  C  (17) 

0  ° 

The  constant  of  integration,  C,  is  eliminated  by  taking  the  vertical 

average  of  equation  (17)  and  subtracting  which  gives  the  following 

expression  for  $: 


C      X 

f-(J*H  dC  -  Jh  dc)  +  ?  (18) 


°  0       0 
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In  the  above  equation  and  throughout  the  remainder  of  this  paper 
represents  the  vertical  average  of  a  given  variable  and     represents 

the  horizontal  average. 

i 

Before  a  solution  for  $  can  be  found,  a  value  for  $  must  be  obtained. 
This  is  accomplished  by  first  taking  the  vertical  average  of  equation  (10) 
thereby  eliminating  the  vertical  advection  term.   Next  the  vertically 
averaged  equation  is  differentiated  with  respect  to  x  yielding: 

The  first  term  of  this  equation  is  zero  since  the  vertical  average  of 
equation  (14)  vanishes  and  the  resulting  equation  is: 


a  fm  3i\    a 


c£ 


(^■iCs^f-^s^ 


o 


(20) 


This  equation  can  be  solved  diagnostically  for  $. 

To  obtain  >—  equation  (18)  is  differentiated  with  respect  to  y 
yielding : 

°    0         0 
Since  ■jr—  =  0  for  this  study 

|  ■  f  <"> 

and  this  quantity  is  assumed  to  be  independent  of  x  to  keep  the 
perturbations  independent  of  y. 

Now  consider  the  domain  averages  of  equations  (10)  and  (11)  which 
are  : 


3F^?  =  -zf  +  f-™!f  +  ^  «*» 
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I_.„  I 


d  . 1    -  d$ 


^  . 


Initially  -~r—  is  computed  geostrophica lly  as  follows 


If  this  relation  is  used  for  all  time  and  if  the  mass  flux  in  the  y 
direction  is  zero  initially,  then  M  will  be  zero  for  all  time.   The 

y 

first  two  terms  on  the  right  side  of  equation  (23)  represent  the  mountain 
drag .   The  sum  of  these  two  terms  was  shown  by  computation  to  be  negative 
in  the  mean  as  expected  from  physical  considerations.   This  will  lead  to 
a  continual  loss  of  x-momentum.   A  loss  such  as  this  would  make  a  quasi- 
steady  final  state  impossible.   An  artificial  source  of  momentum  is  thus 
required  in  this  simple  model  to  maintain  the  mean  zonal  flow.   The 
artificial  source  of  momentum  is  provided  by  holding  ^— ,  as  obtained  in 
equation  (25),  constant  for  all  time.   In  this  case  if  M  is  decreased 
below  the  initial  geostrophic  value  M  will  increase  and  restore  the 

y 

momentum  loss.   Thus,  in  the  long  term  mean,  the  mountain  terms  are 
balanced  by  the  mean  Coriolis  force  due  to  a  net  y-momentum  flux. 
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III.   FINITE  DIFFERENCE  EQUATIONS 

Equations  (10)  through  (14),  and  (18)  are  solved  by  introducing 
finite  differences  in  x,  Q ,   and  t.   The  space  finite  differences  employ 
a  variation  of  a  scheme  developed  by  Arakawa  [1966]  (see  also  Lilly 
[1965]).   The  basic  scheme  was  designed  to  conserve  total  energy  in 
addition  to  conserving  the  squares  of  u,  v,  and  9  with  respect  to  the 
advection  terms.   This  scheme  was  also  designed  to  avoid  the  non- linear 
computational  instability  which  was  discovered  by  Phillips  [1959]. 

The  time  differencing  scheme  employs  a  variation  of  the  Matsuno 
[1966]  two  stage  differencing  scheme  which  is  also  referred  to  as  the 
Euler  simulated  backward  time  step  differencing  scheme.   This  scheme  is 
similar  to  that  used  in  the  general  circulation  model  developed  by  A. 
Arakawa  (see  Langlois  and  Kwok  [1968]).   The  modification  makes  it 
possible  to  avoid  a  checkerboard  spatial  pattern  by  alternating  between 
centered  and  uncentered  space  differences  in  certain  terms. 

The  Euler  backward  difference  scheme  may  be  written  in  the  following 
form 

s*  -  st  +  AtFfc  (step  1) 

st+1  =  sfc  +  ATF*         (step  2) 

where 

s  =  u,  v,  0  and 

F  is  the  right  hand  side  minus  the  linear 
operator  £ (u)  as  defined  by  equations  (10) 
through  (12). 

Centered  space  differences  are  employed  to  compute  F   in  step  1  above 

and  uncentered  space  differences,  using  the  values  of  s*   obtained  in 

step  1,  are  used  to  compute  F*  in  step  2.   During  computation,  forward 
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spatial  differences  are  used  for  even  time  steps  and  backward  spatial 
differences  are  used  for  odd  time  steps.   This  method  allows  geostrophic 
adjustment  to  occur  at  the  smallest  scale  recognized  by  the  model. 

The  vertical  domain  is  divided  into  18  layers  such  that  A£  =  1/18 
for  each  layer.   The  boundary  surfaces  are  denoted  by  the  subscript  i, 
where  i  =  1  at  the  lower  boundary,  the  surface,  to  i  =  19  at  the  top  of 
the  model.   The  quantities  M  ,  M  ,  H,  and  $  are  computed  on  the  even 
levels  and  Q    is  computed  on  the  odd  levels.   Throughout  the  remainder 
of  this  study  j  denotes  the  horizontal  grid  index  with  Ax  as  the  grid 
mesh  length.   For  the  remainder  of  this  paper  CDF,  FDF  and  BDF  will 
represent  centered  difference  form,  forward  difference  form  and  back- 
ward difference  form,  respectively. 

The  finite  difference  form  of  the  linear  operator   (s),  where  is  is 
any  dependent  variable,  can  be  written 

*.;      .  Km  +  si.i><Mx+>-<si.i  +  si.i-i><Mx->]  + 

^S;i,j  2Ax 


In  equation    (26)    the    flux   forms  M       and  M        are   defined  as 

(M      ,    ...  +  M      .     .) 

«      -      x  i,j+l        x  i,y 

M    ,    =  2-J — ^ *-*—     and 

x+  2 

(M      .    .  +  M      .     .    ,) 

X     1,1  X     1,1-1         ,„^„s 

M       ■  *-■-= lJ (CDF) 

x-  2 

M   .    =  M      .     .  ,  .      and 
x+  x    1,J+1 

M       -  M      .         (FDF)   and 
x-  x    l,j 

M^    =  M  and 

x+  x    l,j 

Mx-=Mx   i,j-l    <BDF> 
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(26) 


The  finite  difference  approximation  to  the  integral  in  equation 
(17)  is: 

rK        i=k-2 

H  dC  =   S  (H     +  H   )  K   +  C,  k  >  3  (27) 

i=l    L~*-  >  J     -1-  >  J 
0 

The  value  of  the  constant  of  integration,  C,  is  not  important  since  it 

will  vanish  when  the  vertical  average  of  equation  (27)  is  taken. 

When  i  is  odd,  equation  (13)  can  be  approximated  with: 

$._,_.  .  =  $.  ,  .  +  2AC  (f->)  H.  .  (28) 

1+1, J    i-I, j    ^  \9Qy     i,j 

The  continuity  equation,  equation  (14),  can  be  approximated  with: 

The  flux  forms  M  ,  and  M   are  as  defined  above. 
x+      x- 

Since  all  of  the  dependent  variables  occur  at  each  grid  point,  the 
x  derivatives  are  approximated  by: 

(M)    =  ¥■  (30) 

Vdxy.  .Ax 

The    form  of  the   dependent  variable      6s    is   defined  as  : 

(s .    ...    -   s.    .    .) 
6s   =        L'J+1  2      ^J'1      (CDF) 

6S   "   <Sl,j+l   -   8l,j>    (FDF)   and 
6s   -   (s.    .    -   s.         .)    (BDF) 

To  maintain  consistency  within   the  model   the    finite  difference    form 
of  equation    (20)    is    obtained    from  the    finite   difference    form  of  equation 
(12)   and  appears    in   the    following   form: 
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where 


A4  $   .  -  B.  $.  +  C.  $.  .  -  -  ~  (31) 

J   J+l    J   J    J  J-l     Ax  ' 


B.  =  A.  +  C. 
J    J    J 

A.  ■  Z.,,  and  C.  =  Z.  ,  (used  with  CDF) 
J    J+l      J    J-l 

A.  =  Z.  and  C.  =  Z.  .  (used  with  FDF) 

J    J      J  J-l 

A.  =  Z...  and  C.  =  Z.  (used  with  BDF)  and 
J    J+l      J    J 


i    r 


i  s&s     -j_i  -r  u.  .        U.  .  "T   U.  .  -     >. 


0   i,j   Ax     y  i,j 
o 


where : 


(D  -  -  D   ) 

6D  -  3i=~ — 1— —  (used  with  CDF) 

6D  =  D.  -  D.  ,  (used  with  FDF),  and 
J    J-l 

6D  =  D..,  -  D.  (used  with  BDF) 
J+l    J 

In  the  above  equations  M  , ,  M   ,  and  6s  are  as  defined  above. 

A  solution  for  $  can  now  be  obtained  using  the  exact  technique 
described  on  page  198  of  Richtmyer  and  Morton  [1967]. 
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IV.      SCALE  ANALYSIS 

The   following   scale   analysis    is   done    in   the   z-coordinate   system  as 
a  matter   of  convenience.      The   analysis  will  be   useful    in  determining 
the    initial  conditions    for   this  model  and  will   serve   as   an  aid    in   the 
interpretation   of  the   numerical  solutions. 

The   following  non-dimensional   independent   variables   are  defined: 


x'   =  —  and 


z'  =4 


(32) 


d  (33) 

L  and  d  are  the  horizontal  and  vertical  scales  of  the  disturbance. 
They  are  due  to  the  characteristics  of  the  mountain  and  other  physical 
parameters . 

u  =  U  +  % >  Uu '  '  ,,, . 

o  (34) 

V  "  So  "v'  (35) 

W  =  UB  f  (36)  ' 

do 

9  ■  es  +  B  ar  e"  (37) 

°  o 

The  constant,  U,  represents  the  zonal  flow  away  from  the  mountain.   The 
quantity  ^,  represents  the  stratification  of  the  basic  undisturbed 
atmosphere  which,  for  this  study,  is  assumed  to  be  constant.   The  scale 
of  the  meridional  component  of  the  wind  is: 

V  =  H2 

d  (39) 
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Two  dimensionless  numbers  appear  in  the  analysis.   They  are: 

Ro  "  fL  and  (40) 

f  .£.  .  B 

*o        fL  ~  d  (41) 

R0  is  the  Rossby  number  based  on  the  zonal  wind  and  the  mountain  scale. 
£0  is  the  disturbance  Rossby  number  and  is  so  called  because  it  is  based 
on  the  amplitude  of  the  v  component  of  the  wind. 
The  scale  analysis  requires  that: 


o 
as  long  as  d  <  D. 


(42) 


Using  the  above  definitions  the  steady  state  non-dimensional 
equations  may  be  written: 

£  +  -■+,*(.■#  +  ..  $).„ 


cXi1   bw' 


a^  =  e' 


(44) 
(45) 
(46) 
(47) 


The   boundary  conditions    to  be  applied  are: 
W     ~  ■cbTr  at   2   =  z      and 


(48) 


w'   =0atz'=^ 

d  (49) 
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This  scaling  is  Cnno. 

S  is  consistent  with  th 

—  —.  Note  that  u  th;;rrlcai  -«~ — wm  be 

—  ostrophic   equaaons    ^^         o     -dSo>  e   £Uminated>    ^ 

y  Cnarney   fl962J   ^^ 
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V.   INITIAL  CONDITIONS 

The  profile  of  the  mountain  used  in  this  study  is  given  by  the 
cosine  squared  of  the  horizontal  distance.   The  mountain  is  1200  km 
wide,  has  an  amplitude  designated  B  and  is  duplicated  every  6000  km. 
The  lower  boundary  of  the  model  is  flat  and  horizontal  except  where  the 
mountain  exists. 

The  initial  state  is  derived  from  the  steady  state,  quasi-geostrophic 
equations  which  were  shown  to  be  linear  in  the  preceding  section.   These 
equations  are: 

V  =  |L  (50) 

.•--gr  (5D 

i£ +  w'  ■  °  <52) 

&  +  W-°  <53) 

gr-Q-  (54) 

The  equation  necessary  to  obtain  the  initial  fields  of  u,  v,  and 
9  is  the  potential  vorticity  equation  for  this  set  of  linearized 
equations.   To  obtain  the  potential  vorticity  equation,  differentiate 
equations  (51)  and  (52)  with  respect  to  x.   Subtract  equation  (53) 
from  the  differentiated  form  of  equation  (52).   The  resulting  equation 
is  : 


^<Sr>£- 
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Now  differentiate  equations  (51)  and  (52)  with  respect  to  x  and  equation 
(54)  with  respect  to  z  and  substitute  into  equation  (55).  The  resulting 
equation  is  the  potential  vorticity  equation: 

ax-  C  ,2  +  .  ,2^-  °  (56) 

dz     ox 
To  obtain  the  initial  solution  for  this  study  the  following  boundary 
conditions  are  applied: 

|~-  =  0  at  z'  =  ?  and  (57) 

dz  a 

g;  =  zs  at  z  =  0  (58) 

Equation  (56)  is  satisfied  if: 

^  +  ^  =  0  (58) 

A  general  solution  of  equation  (58),  here  x'  =  kx,  is: 

$'  (x',zf)  =  F(x^)  cos  x'  (59) 

Applying  this  solution  to  equation  (58)  at  the  upper  boundary  and 
solving  for  F(z*)  yields: 

F(z')  =  A  cosh  (z*  -  £)  (60) 

From  this  it  follows  that: 

§'  (x',z')  =  A  cos  x'  cosh  (z*    -  £)  (61) 

The  mountain  topography  can  be  Fourier  analyzed  into  a  cosine 
series.   Each  component  of  that  series  will  have  the  following  form 
if  x  is  non-dimensionalized  with  the  wave  number  of  the  component. 

(62) 
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When  the  solution  for  this  component  is  obtained,  the  quantities  can 

be  dimens ionalized  and  summed  to  satisfy  the  complete  boundary  condition, 

{Differentiation  of  equation  (61)  with  respect  to  z'  and  evaluation 
of  the  resulting  equation  at  z '  =  z'   gives  the  following  result: 

(Ml)      -Acosx'  sinh(z's  -j)  (63) 

z  =z 

s 

If  it  is  assumed  for  this  analysis  that  z'   is  sufficiently  small  that 

equation  (63)  can  be  applied  at  z'   =0,  then  this  equation  can  be 

s 

solved  for  A: 

A  =  -  — i«  (64) 

s inh  I  — 


The  final  solution  for  §'  is 


f  •  =  — ^r   cos  x'  cosh  (z«  -   £)  (65) 

sinh  (  t 


From  this  solution  for  $'  follow  the  solutions  for  u',  v',  and  9': 

1 


^Dx  sin  x'  c 
sinh  (  — 


,  i_ 

u  = p-^r  ' 

sinh  (  t 


osh  (b»  -  £)  (66) 

h  (*»-£)  (67) 


0'  =  ygr  cos  x'  sinh  (z  '  -  £)  (68) 

sinh  (  — 


A  Fourier  series  of  the  full  lower  surface  yields: 

z   =  S  B(k)  cos  kx 
S   k 

in  dimensional  form.   Solutions  (66)  through  (68),  when  dimens iona lized , 

can  be  summed  to  give  the  following  initial  conditions: 
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,       3e 
k  UB<k>  e~  3T 


nw           U  dz 
s  2 

k=i  f2  a  sinh  (aD) 


u  =  U  +  S  ~2 ° cosh  [a(z-D)]  cos  (kx)     (69) 


d0 

V  =  =  f  sinh  (So)   C°Sh  [a(2"D)]  Sin  (kx)  (70) 

k=l 


a9s 
nw  B(k)  ^-^ 

0  =  0  +  E   .  ,  ,*   sinh  [a(z-D)]  cos  (kx) ,  and     (71) 
s       sinh  (aD) 


B   .  /"nw  npm  rf\ 
—  sin  I c ) 

B(k)  =  2 V   "P   <  (72) 


t1"^)] 


where 


nw 

np 


ae    1/2 


«-i(S-sr)     and 


nw  is  the  number  of  wave  numbers  used,  npm  is  the  number  of  points  in 
the  mountain  and  np  is  the  total  number  of  horizontal  grid  points  in 
the  model. 

The  first  20  wave  numbers  were  used  to  initialize  the  wind  and 
temperature  fields.   Any  number  up  to  and  including  20  could  have  been 
used  with  k  =  20  yielding  the  smoothest  initial  solution.   Any  number 
greater  than  20  would  have  introduced  aliasing  since  the  waves  generated 

would  have  a  wave  length  less  than  two  grid  lengths. 

-1 
In  the  initial  fields  the  constants  used  were  U  =  8  m  sec   , 

*r~  =  4.0  °K  km"    ,   D  =  9  km,  |-  =  0.0327  m  sec"      (°K)      ,   and   f  =   10 

'    -1  ° 

sec      .      The  value   of  9     was   computed    for   each   point   from  the   relation: 


d9      r 
G     =  _i 

s       3F 


i*.« 


(C-i)] 
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VI.   RESULTS 

In  this  study  Ax  =  150  km,  At  =  20  min,  and  nine  equally  spaced 

d9  . 

s     o    —  1 
levels  in  zeta  were  used.   Seven  values  of  B  with  ^ —  =4  K  km   and 

s 
eleven  values  of  tt —  at  B  =  1575  m  were  employed. 

Values  of  the  mountain  height  were  computed  using  the  function 

2 

B  cos   (kx) .   This  function  provides  a  smooth  profile  at  the  base  of 

the  mountain.   First,  solutions  were  computed  using  this  function  to 
determine  the  mountain  profile.   The  solutions  were  then  computed  using 
a  linearly  smoothed  mountain  profile.   The  smoothing  was  accomplished 
by  applying  the  following  linear  formula: 

V 


+  Z.  +  Z. 
***■      ( -JZL 
'3 


z*-i,p+i  •  j  ■  -j-i)  (73) 


This  procedure  was  repeated  four  times.   The  resulting  mountain  height 
was  about  two  thirds  that  of  the  original  mountain  height  and  the 
horizontal  extent  was  approximately  doubled.   In  each  case  the  solutions 
were  allowed  to  converge  toward  a  quasi-steady  state.   To  permit  com- 
parison of  the  various  cases  the  results  were  converted  from  the  zeta- 

surfaces  to  the  z-surfaces  at  one  km  intervals.   In  this  study  R   is 

o 

independent  of  the  mountain  height  and  its  value  is  R  =  .419. 

The  general  features  of  the  flow  were  basically  the  same  for  all 
cases  and  were  similar  to  those  reported  by  Walton  [1968].   There  was 
a  ridge  in  the  pressure  field  over  the  mountain.   The  v  component  of  the 
wind  was  from  the  south  on  the  windward  side  of  the  mountain  and  from 
the  north  on  the  leeward  side.   The  u  component  of  the  wind  was  at  a 
maximum  near  the  surface  over  the  mountain  and  decreased  rapidly  with 
height  to  a  minimum  at  the  top  of  the  atmosphere.   Over  the  plain, 
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away  from  the  mountain,  the  u  component  was  a  minimum  near  the  surface 
and  increased  very  slightly  with  height  to  a  maximum  at  the  top  of  the 

atmosphere.   The  potential  temperature  was  a  minimum  over  the  mountain 

! 

crest . 

Fig.  2  is  a  graph  of  the  difference  between  the  maximum  value  of  $ 
in  the  ridge  and  the  minimum  value  of  $  at  the  bottom  of  the  trough,  at 
the  three  km  level,  with  respect  to  mountain  height.   The  upper  curve 
represents  the  smoothed  solutions  from  the  unsmoothed  mountain. 

The  difference  between  the  smoothed  and  unsmoothed  solutions  is 
least  for  mountain  heights  less  than  900  m.   Below  this  height  the 

solutions  are  almost  identical.   The  scale  analysis  of  the  basic 

2 

equations  used  in  this  study  showed  that  when  R   is  small  the  balance 

approximation  should  hold.   When  Si      is  also  small  the  equations  reduce 
to  the  linear  quasi-geostrophic  equations.   Since  the  equations  are 
linear  f  or  $   small  the  solutions  should  be  nearly  equal  in  this  range 
as  they  actually  are.   At  mountain  heights  greater  than  900  m  the 
non-linear  effects  become  large  enough  to  cause  a  marked  separation 
of  the  two  curves.   These  results  are  comparable  in  magnitude  to  those 
obtained  by  Walton  [1968].   In  his  study  he  applied  a  different  differ- 
encing scheme  to  obtain  solutions  to  this  same  model. 

Theaddition  of  the  artificial  source  of  momentum  suggested  in  Section 
II  appears  to  have  been  valid.   The  slow  decay  in  mean  zonal  velocity 
with  time  noted  by  Walton  [1968]  did  not  occur.   In  providing  the 
artificial  source  of  momentum  it  was  assumed  that  the  model  would  come 
into  adjustment  and  that  the  mountain  drag  term  would  be  balanced  by 
the  term  denoting  mass  flux  in  the  y-direction.   Fig.  3  shows  a  plot  of 
these  two  quantities  versus  time.   This  figure  shows  the  values  for  a 
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single  mountain  height,  however  it  is  representative  of  all  the  mountain 
heights  in  this  study.   The  upper  curve  represents  the  momentum  flux  in 

the jy-direction  and  the  lower  curve  represents  the  mountain  drag.   It 

i 

can  be  seen  that  these  two  quantities  do  behave  as  expected  after  the 
initial  period  when  the  model  is  coming  into  adjustment.   It  should  be 
noted  that  the  sum  of  the  two  curves  does  not  equal  zero.   The  long 
term  change  in  the  x-mass  flux  is  very  small.   The  reason  for  this 
difference  is  apparently  due  to  the  fact  that  centered  space  differences 
were  used  to  calculate  the  mountain  drag  while  uncentered  differences 
were  used  to  calculate  the  actual  forecast  variables. 

The  solutions,  represented  by  A$ ,  did  not  appear  to  be  converging 
significantly  at  the  end  of  40  days  into  the  forecast  period.   For  the 
higher  values  of  mountain  height  there  appeared  to  be  slight  conver- 
gence, however  this  may  have  been  due  to  the  time  finite  differencing 
scheme  which  tends  to  damp  the  computational  mode  of  solution.   Fig.  4 
is  an  example  of  the  oscillation  of  the  solution  obtained.   These' 
oscillations  have  a  period  of  about  eight  days.   This  pattern  was 
identical  for  all  cases  with  only  the  amplitude  of  oscillation  varying. 
The  amplitude  increased  with  increasing  mountain  height.   This  oscillation 
appears  to  be  the  result  of  mechanical  interaction  between  the  mountains. 
One  case  involving  a  mountain  height  of  2025  m  was  computed  retaining 
the  width  of  the  mountain  at  1200  km  but  doubling  the  width  of  the  plain 
from  4800  km  to  9600  km.   As  a  result  of  this  the  number  of  oscillations 
occurring  in  a  40  day  forecast  was  diminished  by  one  half  with  the  ampli- 
tude of  oscillation  decreasing  slightly.   The  values  used  in  Fig.  2  are 
the  mean  values,  at  a  given  height,  about  which  A$  oscillated. 

As  a  check  of  the  validity  of  the  geostrophic  assumption  the  v  com- 
ponent of  the  geostrophic  wind  was  computed.   Except  in  the  near  vicinity 
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of  the  mountain  the  values  obtained  from  the  forecast  equations  and  the 
values  obtained  from  the  geostrophic  computations  were  similar.   The 
difference  between  the  geostrophic  and  actual  v  component  was  less  than 
ten  percent  in  conformity  with  equation  (13). 

The  variation  of  mountain  drag  with  mountain  height  and  with 
static  stability  was  also  investigated.   Fig.  5  displays  the  results 
of  the  variation  of  mountain  drag  with  mountain  height.   From  this 
figure  it  can  be  seen  that  mountain  drag  increases  quadratically  with 
increasing  mountain  height.   Fig.  6  shows  the  relationship  between 
mountain  drag  and  static  stability.   This  figure  shows  that  mountain 
drag  increases  linearly  with  increasing  static  stability.   An  investi- 
gation of  the  gross  features  of  the  flow  for  increasing  static  stability 
at  a  constant  mountain  height  showed  no  significant  changes  except  in 
magnitude . 

The  stream  function  was  computed  and  streamlines  were  drawn  for 
each  mountain  height  at  selected  time  spacing.   Fig.  7  is  a  typical 
example  of  the  streamlines  obtained.   It  shows  horizontal  flow  over 
the  plain  at  all  levels  and  flow  conforming  to  the  mountain  in  the 
vicinity  of  the  mountain  with  the  upper  level  flow  nearly  horizontal 
everywhere . 

In  the  case  where  B  =  2250  m  static  instability  developed  about 
36  hours  into  the  forecast  period.   A  comparison  of  the  u  and  $  fields 
obtained  in  the  hydraulic  jump  case  B  described  by  Houghton  and 
Kasahara  [1968]  with  those  obtained  in  this  study  for  B  =  2250  m  shows 
a  marked  similarity  in  the  structure  of  the  fields.   Houghton  and 
Kasahara  used  a  single  layer,  non-rotating  model.   The  lower  values  of 
mountain  height  used  in  the  present  study  apparently  fall  into  the 
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non-critical,  or  no  jump  domain,  described  by  Houghton  and  Kasahara. 
As  mountain  height  increases  the  values  eventually  fall  in  the  hydraulic 
jump  domain.   The  reason  for  the  static  instability  at  B  =  2250  m  is 
thus  attributed  to  hydraulic  jump.   Fig.  8  shows  the  values  of  u 
initially  and  six  hours  later  for  B  =  2250  m  at  the  three  km  level. 
The  development  of  the  minimum  in  u  in  the  lee  of  the  mountain  is 
already  evident  even  though  static  instability  has  not  yet  occurred. 
This  model  is  not  appropriate  for  the  treatment  of  hydraulic  jump 
because  of  approximations  which  have  been  made.   Calculations  of 
mountain  drag  were  not  made  for  this  case  but  it  is  expected  that  if 
these  calculations  had  been  made  the  character  of  the  drag  properties 
would  have  been  modified  above  B  =  2025  m. 
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VIi.  conclusions 


h  wed  that  linear  smoothing  of  the  topography 
obtained  in  this  study  showed  that 

=«»,  as  the  disturbance  Rossby  Number 
an  error  which  Increases  as 

,-   Se„^   1  „  fh„  results  from  this  study 

r     B   B_  (% SA     increases.   If  the  resuiu 

|A,  =  d  =  f  L  V.9   °z7   J  R,  ls  larger ,  the 

could  be  extrapolated  to  smaller  scales,  I.e.,       0 
dl££erence  between  results  obtained  using  a  smoothed  topography 
those  obtained  using  an  nnsmoothed  topography  would  be  even  greate, 
L  is  consistent  with  the  example  given  by  Charney  [lW,  U  -h 
he  considered  an  inf initesimally  thin  mountain,  oriented  along  a 

k0T-0   Without  smoothing 
a-    „   hn  the  top  of  the  atmosphere,   wuno 
meridian  and  extending  to  the  top 

k   lute  barrier,  however,  with  smoothing  its 
it  would  act  as  an  absolute  barrie  , 
vertical  extent  would  be  decreased  to  nothing  and  it  would  have  no 

•   •   «„«-  rlear  what  error  will 
hpric  flow.   However  it  is  not  clear 
effect  on  atmospheric  now. 
be  found  if  becomes  large  enough  for  pressure  .umps  to  £otm  ^ 

error  behavior  will  probably  be  different  than  that  obtained  here 

^  addition  of  an  artificial  source  of  momentum  did  maintain  the 

onent  of  the  flow  at  a  value  near  the  initral 
-n  value  of  the  u  component  ^ 

.,  zonal  flow.  A  comparison  of  the  values    *, 

nofiftl  did  not  show  any  significant 
,   u    „ufflined  bv  Walton  [19bbJ  aia 
study  and  those  obtained.  Dy  « 

,       ta  the  results  obtained  for  smooth  and  nnsmoothed  mountams . 
difference  In  the  re&uj-L. 

f  the  olain  between  the  mountains  it  was 
By  doubling  the  size  of  the  plain 

,   •    f   t-hP  M>    field  was  reduced  by 
observed  that  the  period  of  oscillatton  of  the 

i  A„A   that  this  oscillation  is  due,  at  least 
one  half.   It  was  thus  concluded  that 
U  Part,  to  mechanical  interactions  between  the  mountains. 


Ul 


Mountain  drag  computed  for  varying  mountain  heights  at  a  constant 
static  stability  showed  an  apparent  quadratic  relationship  between 
increasing  drag  and  increasing  mountain  height.   Mountain  drag  computed 
for  varying  static  stability  at  a  constant  mountain  height  showed  an 
apparent  linear  relationship  between  increasing  drag  and  increasing 
stability. 

For  the  case  when  B  =  2250  m  static  instability  resulted  early  in 
the  forecast  period.   The  reason  for  this  instability  was  attributed 
to  the  formation  of  a  hydraulic  jump.   This  cannot  be  definitely  stated 
as  the  model,  in  the  form  which  was  used  for  this  study,  was  not 
appropriate  for  the  treatment  of  hydraulic  jump  because  of  the  lack  of 
friction  and  because  of  certain  approximations  which  had  been  made. 

Further  studies  with  this  model  should  be  conducted  using  a  wider 
range  of  mountain  scales.   This  would  give  a  more  complete  relation- 
ship between  mountain  drag  and  mountain  structure  which  could  be  used 
in  parameterizing  mountain  effects  in  numerical  weather  prediction 
models.   The  drag  could  be  obtained  for  the  larger  values  of  $  if 
appropriate  treatment  of  the  hydraulic  jump  is  included.   The  mechanical 
interaction  between  the  periodic  mountains  could  be  reduced  if  a  variable 
mesh  length  grid  in  the  horizontal  scale  of  the  model  is  used.   This 
would  permit  a  short  mesh  length  in  the  vicinity  of  the  mountain  where 
most  of  the  activity  is. 
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